CCD Project Results created by Jordan Muraskin

In [43]:
% pylab inline
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot,plot
import matplotlib.pyplot as plt
init_notebook_mode()
import CCD_packages
reload(CCD_packages)

Populating the interactive namespace from numpy and matplotlib
Out[43]:
<module 'CCD_packages' from 'CCD_packages.py'>
In [42]:
! git pull

remote: Counting objects: 4, done.
remote: Compressing objects: 100% (1/1), done.
remote: Total 4 (delta 3), reused 4 (delta 3), pack-reused 0
Unpacking objects: 100% (4/4), done.
From https://github.com/jordanmuraskin/CCD-scripts
   2d52ed2..4506719  master     -> origin/master
Updating 2d52ed2..4506719
Fast-forward
 analysis/CCD_packages.py | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

Get Subject Info

In [3]:
GroupDF,motionInfo=CCD_packages.getCCDSubjectData()

Select which subjects to use

In [4]:
goodsubj = CCD_packages.getSubjectList(GroupDF=GroupDF)

Resting State Network Time Series with Model

In [31]:
CCD_packages.createTimeSeriesPlots(GroupDF,goodsubj,DMN_name='RSN3',title='Default Mode Activity',ylabel='DMN Activity',figsize=(9,4.5))
CCD_packages.createTimeSeriesPlots(GroupDF,goodsubj,DMN_name='RSN0',title='Visual IC Activity',ylabel='Visual Activity',figsize=(9,4.5))

Model Correlations for Each Subject

In [6]:
CCD_packages.createSubjectModelBarPlot(GroupDF,goodsubj,figsize=(9,4.5),savefig=True)

Scan Order model correlations

In [38]:
CCD_packages.createScanOrderBarPlot(GroupDF,goodsubj)

Output Mean Subject Model Correlations for Default Mode Network

In [8]:
# print 'RSN3'
CCD_packages.printModelCorrelations(GroupDF,goodsubj,'RSN3')

No Feedback Focus Correlation= -0.70
Feedback Focus Correlation= -0.57
No Feedback Wander Correlation= 0.67
Feedback Wander Correlation= 0.59
No Feedback Overall Correlation= 0.70
Feedback Overall Correlation= 0.60
In [9]:
print 'RSN0'
CCD_packages.printModelCorrelations(GroupDF,goodsubj,'RSN0')

RSN0
No Feedback Focus Correlation= -0.54
Feedback Focus Correlation= -0.67
No Feedback Wander Correlation= 0.55
Feedback Wander Correlation= 0.66
No Feedback Overall Correlation= 0.56
Feedback Overall Correlation= 0.69
In [10]:
hmFB,hmNFB,hmDiff=CCD_packages.generateHeatMaps(GroupDF,goodsubj)

Resting State Network Connectivity Matrices

Taken from Dual Regression Outputs with PNAS RSN10

In [13]:
fig=CCD_packages.heatmap2Chord(hmFB,plotName='FeedbackChordDiagram',title='Network Correlations with Feedback On',savefig=False,scale=[-1,1])
iplot(fig)
In [14]:
iplot(CCD_packages.heatmap2Chord(hmNFB,plotName='FeedbackChordDiagram',title='Network Correlations with Feedback Off',savefig=False,scale=[-1,1]))
In [15]:
from scipy.stats import ttest_1samp
from mne.stats.multi_comp import fdr_correction
from numpy import triu_indices,tril_indices
t,p=ttest_1samp(hmDiff,0,axis=2)
row,column=triu_indices(10,1)
rowl,columnl=tril_indices(10)
# sns.heatmap(t)
p05,padj=fdr_correction(p[row,column],0.05)
fdr_corrected=t.copy()
fdr_corrected[rowl,columnl]=0
fdr_corrected[row,column]=t[row,column]*p05

Pairwise Differences between Feedback On vs. Feedback Off Resting State Network Connectivity

In [17]:
iplot(CCD_packages.heatmap2Chord(fdr_corrected,plotName='FeedbackDiffChordDiagram',title='Network Correlations differences between Feedback On/Off',savefig=False,scale=[-10,10]))
In [19]:
predictions,coefs,performance,fb_coefs,nfb_coefs=CCD_packages.linearRegressionData(GroupDF,goodsubj,numFolds=10)

Running Feedback on Regressions
Finished...
Running Feedback off Regressions
Finished...

10-Fold Linear Regression Model Prediction with RSN Timeseries

In [39]:
CCD_packages.createRegressionPlots(predictions,performance,coefs,fb_coefs,nfb_coefs,GroupDF,goodsubj)
  • *-Indicates p<0.05 FDR-Corrected Feedback On vs. Feedback Off
In [30]:
f, axarr = plt.subplots(2, sharex=True,figsize=(16,8))

TFCEposImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/Feedback/motionThresh-0.200000/cope1/cope1/cope1_tfce_corrp_tstat1.nii.gz'
posImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/Feedback/motionThresh-0.200000/cope1/cope1/cope1_tstat1.nii.gz'
TFCEnegImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/Feedback/motionThresh-0.200000/cope1/cope1/cope1_tfce_corrp_tstat2.nii.gz'
negImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/Feedback/motionThresh-0.200000/cope1/cope1/cope1_tstat2.nii.gz'
CCD_packages.createTFCEfMRIOverlayImages(TFCEposImg,posImg,TFCEnegImg,negImg,title='Focus versus Wander with Feedback On',plotToAxis=True,axes=axarr[0],f=f)

TFCEposImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/noFeedback/motionThresh-0.200000/cope1/cope1/cope1_tfce_corrp_tstat1.nii.gz'
posImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/noFeedback/motionThresh-0.200000/cope1/cope1/cope1_tstat1.nii.gz'
TFCEnegImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/noFeedback/motionThresh-0.200000/cope1/cope1/cope1_tfce_corrp_tstat2.nii.gz'
negImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/noFeedback/motionThresh-0.200000/cope1/cope1/cope1_tstat2.nii.gz'
CCD_packages.createTFCEfMRIOverlayImages(TFCEposImg,posImg,TFCEnegImg,negImg,title='Focus versus Wander with Feedback Off',plotToAxis=True,axes=axarr[1],f=f)

# f.tight_layout()
Out[30]:
<nilearn.plotting.displays.ZSlicer at 0x7fbafc22c490>
In [29]:
f, axarr = plt.subplots(1,figsize=(16,4))
TFCEposImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope1/cope1/cope1_tfce_corrp_tstat1.nii.gz'
posImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope1/cope1/cope1_tstat1.nii.gz'
TFCEnegImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope1/cope1/cope1_tfce_corrp_tstat2.nii.gz'
negImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope1/cope1/cope1_tstat2.nii.gz'
CCD_packages.createTFCEfMRIOverlayImages(TFCEposImg,posImg,TFCEnegImg,negImg,vmax=5,title='Feedback versus No Feedback F>W',plotToAxis=True,f=f,axes=axarr)
# f.tight_layout()
Out[29]:
<nilearn.plotting.displays.ZSlicer at 0x7fba6848ebd0>
In [28]:
f, axarr = plt.subplots(3, sharex=True,figsize=(16,12))

names=['Focus','Wander','Mean']

for indx,cope in enumerate(range(3,6)):
    TFCEposImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope%d/cope%d/cope%d_tfce_corrp_tstat1.nii.gz' % (cope,cope,cope)
    posImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope%d/cope%d/cope%d_tstat1.nii.gz' % (cope,cope,cope)
    TFCEnegImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope%d/cope%d/cope%d_tfce_corrp_tstat2.nii.gz' % (cope,cope,cope)
    negImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/randomise/paired-Ttest/motionThresh-0.200000/cope%d/cope%d/cope%d_tstat2.nii.gz' % (cope,cope,cope)
    CCD_packages.createTFCEfMRIOverlayImages(TFCEposImg,posImg,TFCEnegImg,negImg,title='%s Feedback vs. No Feedback' % names[indx],plotToAxis=True,axes=axarr[indx],f=f)
# f.tight_layout()
In [32]:
f, axarr = plt.subplots(1,figsize=(16,5))
TFCEposImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/PPI_DMN/randomise/Feedback/motionThresh-0.200000/cope3/cope3/cope3_tfce_corrp_tstat1.nii.gz'
posImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/PPI_DMN/randomise/Feedback/motionThresh-0.200000/cope3/cope3/cope3_tstat1.nii.gz'
TFCEnegImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/PPI_DMN/randomise/Feedback/motionThresh-0.200000/cope3/cope3/cope3_tfce_corrp_tstat2.nii.gz'
negImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/PPI_DMN/randomise/Feedback/motionThresh-0.200000/cope3/cope3/cope3_tstat2.nii.gz'
CCD_packages.createTFCEfMRIOverlayImages(TFCEposImg,posImg,TFCEnegImg,negImg,slices=[-10,41,43,45],vmax=5,title='Feedback DMN PPI',plotToAxis=True,f=f,axes=axarr)
# f.tight_layout()
Out[32]:
<nilearn.plotting.displays.ZSlicer at 0x7fba68635210>
In [34]:
f, axarr = plt.subplots(1,figsize=(16,5))
TFCEposImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/DMN/paired/motionThresh-0.200000/DMN_pair/paired_tfce_corrp_tstat1.nii.gz'
posImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/DMN/paired/motionThresh-0.200000/DMN_pair/paired_tstat1.nii.gz'
TFCEnegImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/DMN/paired/motionThresh-0.200000/DMN_pair/paired_tfce_corrp_tstat2.nii.gz'
negImg='/home/jmuraskin/Projects/CCD/working_v1/groupAnalysis/DMN/paired/motionThresh-0.200000/DMN_pair/paired_tstat2.nii.gz'
CCD_packages.createTFCEfMRIOverlayImages(TFCEposImg,posImg,TFCEnegImg,negImg,slices=range(-10,50,15),vmax=5,title='Dual Regression DMN Feedback versus No Feedback',plotToAxis=True,f=f,axes=axarr)
# f.tight_layout()
Out[34]:
<nilearn.plotting.displays.ZSlicer at 0x7fba6380fd50>

Default Mode Regulation versus Phenotypic Measures

Very Preliminary

First lets use only total scores as our independent variables to predict DMN correlation with Model

In [51]:
from pandas import read_csv
#phenotypic info
phenoFile='/home/jmuraskin/Projects/CCD/Pheno/narsad+vt_new.csv'

pheno=read_csv(phenoFile)
pheno=pheno.set_index('participant').fillna(value=0)
subject_list=goodsubj

addmotion=True
# phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDAIM_42','V1_CCDAIM_43','V1_CCDAIM_44','V1_CCDAIM_45','V1_CCDSIPI_46','V1_CCDSIPI_47','V1_CCDSIPI_48','V1_CCDRSQ_72','V1_CCDRSQ_73','V1_CCDRSQ_74','V1_CCDRSQ_75','V1_CCDPANAS_21','V1_CCDPANAS_22','V1_CCDERQ_11','V1_CCDERQ_12']
phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDAIM_41','V1_CCDSIPI_49','V1_CCDPANAS_23','V1_CCDERQ_13','V1_DBDI_22']
# phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDPANAS_21','V1_CCDPANAS_22']
# phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDAIM_42','V1_CCDAIM_43','V1_CCDAIM_44','V1_CCDAIM_45']
# phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDRSQ_72','V1_CCDRSQ_73','V1_CCDRSQ_74','V1_CCDRSQ_75']
# phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDSIPI_46','V1_CCDSIPI_47','V1_CCDSIPI_48']

modelInfo=np.arctanh(GroupDF[np.all([GroupDF.Subject_ID.isin(subject_list),GroupDF.FB=='FEEDBACK'],axis=0)].groupby('Subject_ID').mean()['modelcorr'])
motionInfo=GroupDF[np.all([GroupDF.Subject_ID.isin(subject_list),GroupDF.FB=='FEEDBACK'],axis=0)].groupby('Subject_ID').mean()['meanFD']

if addmotion:
    pheno=pheno.loc[subject_list][phenoValues]
    pheno['meanFD']=motionInfo
    modelX=pheno
else:
    modelX=pheno.loc[subject_list][phenoValues]
    
inputs=modelInfo
# inputs=array(fb_performance['R'])
# inputs=array(fb_coefs[fb_coefs.pe==3]['Coef'])
results=CCD_packages.runRLMR(inputs,modelX,phenoValues + ['meanFD'],RLM=True)

                    Robust linear Model Regression Results                    
==============================================================================
Dep. Variable:              modelcorr   No. Observations:                   37
Model:                            RLM   Df Residuals:                       28
Method:                          IRLS   Df Model:                            8
Norm:                          HuberT                                         
Scale Est.:                       mad                                         
Cov Type:                          H1                                         
Date:                Thu, 25 Aug 2016                                         
Time:                        16:53:50                                         
No. Iterations:                    46                                         
==================================================================================
                     coef    std err          z      P>|z|      [95.0% Conf. Int.]
----------------------------------------------------------------------------------
const              0.1224      0.216      0.567      0.571        -0.301     0.545
V1_DEM_001        -0.0036      0.002     -2.013      0.044        -0.007 -9.59e-05
V1_DEM_002         0.0151      0.040      0.377      0.706        -0.063     0.093
V1_CCDAIM_41      -0.0013      0.001     -1.481      0.139        -0.003     0.000
V1_CCDSIPI_49      0.0019      0.001      2.115      0.034         0.000     0.004
V1_CCDPANAS_23    -0.0028      0.002     -1.210      0.226        -0.007     0.002
V1_CCDERQ_13       0.0092      0.002      4.099      0.000         0.005     0.014
V1_DBDI_22        -0.0026      0.002     -1.074      0.283        -0.007     0.002
meanFD            -1.4986      0.593     -2.529      0.011        -2.660    -0.337
==================================================================================

If the model instance has been used for another fit with different fit
parameters, then the fit options might not be the correct ones anymore .
In [53]:
phenoValues=['V1_DEM_001','V1_DEM_002','V1_CCDAIM_42','V1_CCDAIM_43','V1_CCDAIM_44','V1_CCDAIM_45','V1_CCDSIPI_46','V1_CCDSIPI_47','V1_CCDSIPI_48','V1_CCDRSQ_72','V1_CCDRSQ_73','V1_CCDRSQ_74','V1_CCDRSQ_75','V1_CCDPANAS_21','V1_CCDPANAS_22','V1_CCDERQ_11','V1_CCDERQ_12']
pheno=read_csv(phenoFile)
pheno=pheno.set_index('participant').fillna(value=0)
subject_list=goodsubj
if addmotion:
    pheno=pheno.loc[subject_list][phenoValues]
    pheno['meanFD']=motionInfo
    modelX=pheno
else:
    modelX=pheno.loc[subject_list][phenoValues]
    
inputs=modelInfo
# inputs=array(fb_performance['R'])
# inputs=array(fb_coefs[fb_coefs.pe==3]['Coef'])
results=CCD_packages.runRLMR(inputs,modelX,phenoValues + ['meanFD'],RLM=False)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:              modelcorr   R-squared:                       0.820
Model:                            OLS   Adj. R-squared:                  0.658
Method:                 Least Squares   F-statistic:                     5.076
Date:                Thu, 25 Aug 2016   Prob (F-statistic):           0.000511
Time:                        16:55:03   Log-Likelihood:                 48.838
No. Observations:                  37   AIC:                            -61.68
Df Residuals:                      19   BIC:                            -32.68
Df Model:                          17                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [95.0% Conf. Int.]
----------------------------------------------------------------------------------
const              0.1508      0.276      0.547      0.591        -0.426     0.728
V1_DEM_001        -0.0054      0.002     -2.337      0.031        -0.010    -0.001
V1_DEM_002         0.0080      0.046      0.176      0.862        -0.088     0.104
V1_CCDAIM_42      -0.0026      0.002     -1.149      0.265        -0.007     0.002
V1_CCDAIM_43       0.0042      0.005      0.814      0.426        -0.007     0.015
V1_CCDAIM_44      -0.0070      0.005     -1.402      0.177        -0.017     0.003
V1_CCDAIM_45       0.0103      0.004      2.598      0.018         0.002     0.019
V1_CCDSIPI_46     -0.0043      0.003     -1.373      0.186        -0.011     0.002
V1_CCDSIPI_47      0.0083      0.002      3.634      0.002         0.004     0.013
V1_CCDSIPI_48     -0.0006      0.002     -0.249      0.806        -0.006     0.004
V1_CCDRSQ_72       0.0020      0.003      0.777      0.447        -0.003     0.007
V1_CCDRSQ_73      -0.0075      0.008     -0.933      0.362        -0.024     0.009
V1_CCDRSQ_74       0.0005      0.007      0.071      0.944        -0.013     0.014
V1_CCDRSQ_75       0.0091      0.007      1.296      0.211        -0.006     0.024
V1_CCDPANAS_21     0.0006      0.004      0.165      0.870        -0.007     0.008
V1_CCDPANAS_22    -0.0137      0.004     -3.330      0.004        -0.022    -0.005
V1_CCDERQ_11       0.0119      0.004      3.355      0.003         0.004     0.019
V1_CCDERQ_12       0.0053      0.005      1.125      0.274        -0.005     0.015
meanFD            -1.8342      0.604     -3.038      0.007        -3.098    -0.570
==============================================================================
Omnibus:                        3.619   Durbin-Watson:                   2.512
Prob(Omnibus):                  0.164   Jarque-Bera (JB):                1.784
Skew:                           0.209   Prob(JB):                        0.410
Kurtosis:                       2.009   Cond. No.                     1.14e+17
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 4.99e-29. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
In [ ]: